//
//  steady.cpp
//  NLCEQ_mixed
//
//  Created by Duong Ngo on 8/6/16.
//  Copyright © 2016 Duong Ngo. All rights reserved.
//

#include <iostream>
#include <vector>
#include <math.h>

#include "steady.hpp"
#include "global.hpp"

void find_steady_state(std::vector<double>& ss)
{
    //Initial Guess
    double cb= 0.00343596;
    double ch= 2.11634218;
    double k= 14.54133172;
    
    //Other variables
    double gamma= 1/cb;
    double mur= gamma*(1-beta_b*Rn);
    double ql= beta_b*(deltab- thetaa)/(1-beta_b*deltab);
    double tau =0 ;
    double i= deltak*k;
    double bh= psi*k;
    double s= (1-deltab)*bh;
    double Rm= (gamma- varphi*mur)/(beta_b*gamma);
    double m = (ch+i+deltab*bh- ql*s)/Rm;
    double Rf= 1/beta_b;
    double n= varphi*m;
    double muc= n+(1-kappa)*bh- m;
    double pm= (epsilon-1)/epsilon;
    double lamb= 1/ch;
    double lama= beta_h*Rm*lamb;
    double etaz= 1/ch- lama;
    double etab= ql*lamb- beta_h*(deltab+(1-deltab)*ql)*lamb;
    double ip=1;
    double pie=1;
    double u= 1/beta_b;
    double y= cb+ ch+ i+ thetaa*bh;
    double l= (1-alphaa)*pm*y*lama/chi;
    
    
    ss[0]=gamma; ss[1]=cb; ss[2]=Rf; ss[3]=Rm; ss[4]=ql;
    ss[5]=mur; ss[6]=muc; ss[7]=ip; ss[8]=pie; ss[9]=n;
    ss[10]=m; ss[11]=bh; ss[12]=k; ss[13]=tau; ss[14]=s;
    ss[15]=ch; ss[16]=etaz; ss[17]=lama;  ss[18]=lamb; ss[19]=etab;
    ss[20]=i; ss[21]=pm; ss[22]=y; ss[23]=u; ss[24]=l;
  
    
}